METHOD AND APPARATUS FOR LOCALIZING 



BIOMAGNETIC SIGNALS 

Technical Field 

5 [0001] The invention relates to the measurement of magnetic fields 
generated by biological activity (biomagnetic signals). The invention has 
application to the localization of sources in magnetic imaging of 
biological structures. Some embodiments of the invention are used to 
provide continuous head localization information in 
1 0 magnetoencephalography . 

Background 

[0002] Magnetoencephalography ("MEG") involves detecting 
magnetic fields produced within a subject's brain. Information about 
1 5 such biomagnetic fields is most useful when it can be associated with 
particular structures within the subject's brain. One can localize the 
subject's head relative to the measured magnetic fields. To do this, the 
position of the subject's head relative to the magnetic detectors used to 
detect the magnetic fields must be known. 

20 

[0003] One way to localize a subject's head is to attach small coils 
at three or more known locations on the subject's head. The coils create 
fluctuating magnetic fields when alternating electrical currents are passed 
through them. Magnetic detectors are used to detect the coils' magnetic 
25 fields. The location of the coils, and thus the location of the subject's 
head, can be determined from the detected magnetic fields of the coils. 
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[0004] One problem with some current head localization systems is 
that the coils generate large magnetic fields that can saturate the very 
sensitive magnetic field detectors used to detect biomagnetic signals. 
Because of this, head localization cannot be performed while 
5 biomagnetic signals are being measured. With such current systems it is 
necessary to perform head localization before and/or after acquiring data 
about the biomagnetic signals. Such systems assume that the position of 
the subject's head changes predictably and can be determined from its 
positions before and/or after the biomagnetic signals are measured. These 
10 assumptions are not always accurate. ^ 

[0005] Acquiring biomagnetic signals is often performed over 
significant time spans. During these times a subject's head is likely to 
move even if the subject is attempting to stay still. 

15 

[0006] Magnetic noise can exacerbate the problems with head 
localization. The magnetic detectors used to detect the magnetic fields 
from the coils will also pick up noise signals, such as signals at the 
power line frequency (generally 60 Hz in North America and 50 Hz in 

20 Europe) and at harmonics of the powerline frequency. Such noise signals 
can degrade the accuracy with which the positions of each of the coils 
can be determined. Magnetic noise includes environmental background 
noise, noise resulting from operation of the magnetic detectors and other 
components in the signals from the magnetic detectors which do not arise 

25 from the magnetic fields of the coils. 
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[0007] There is a need for systems which permit continuous 
localization of heads and other structures during the acquisition of 
biomagnetic signals. There is a particular need for such systems which 
are relatively insensitive to noise. 

5 

Summary of the Invention 

[0008] This invention relates to localizing sources of biomagnetic 
signals. One aspect of the invention provides a method for locating an 
object relative to an array of magnetic sensors in an environment in 

10 which there is present a noise signal having a fundamental frequency 
/noise* the method comprising generating one or more magnetic signals 
by means of one or more magnetic emitters mounted at known locations 
on the object, the one or more magnetic signals having one or more 
frequencies; during an integration time T, detecting the one or more 

15 magnetic signals and the noise signal at six or more magnetic detectors; 
and, determining relative amplitudes of the magnetic signals; wherein 
the one or more frequencies of the magnetic signals are substantially 
equal to frequencies at which a power spectrum of the detected noise 
signal has zeros. 

20 

[0009] Further aspects of the invention and features of specific 
embodiments of the invention are described below. 

Brief Description of the Drawings 
25 [0010] In drawings which illustrate non-limiting embodiments of 
the invention, 



r 
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Figure 1 is a schematic diagram of a MEG system which performs 
head localization according to one embodiment of the invention; and, 

Figure 2 is a flow chart which illustrates a method for locating an 
object relative to an array of magnetic sensors according to another 
5 embodiment of the invention. 



Description 

[0011] Throughout the following description, specific details are 
set forth in order to provide a more thorough understanding of the 
10 invention. However, the invention may be practiced without these 
particulars. In other instances, well known elements have not been 
shown or described in detail to avoid unnecessarily obscuring the 
invention. Accordingly, the specification and drawings are to be 
regarded in an illustrative, rather than a restrictive, sense. 

15 

[0012] Figure 1 shows a MEG system 20 according to one 
embodiment of the invention. System 20 includes a number M of 
magnetic detectors 24 which are located to detect magnetic fields 
originating within the head of a subject S. While Figure 1 represents a 

20 detector array by only a few detectors 24, a MEG system typically has a 
few hundred detectors 24. Detectors 24 may be SQUID detectors. Signals 
from detectors 24 are appropriately conditioned by way of suitable signal 
conditioning circuitry 26. The processed signals are provided to a signal 
processing system 28. Signal processing system 28 generates output at an 

25 output device 30. The output may comprise, for example, a graphical 
display, a data file containing MEG data, or the like. 



[0013] Coils 32, 33 and 34 are attached to the subject's head and 
are respectively driven with signals of frequencies f } ,f 2 and f 3 
(collectively referred to as f g ) by a controller 36 which, in the illustrated 
embodiment includes oscillators 37, 38, and 39. Driving signals for coils 
5 32, 33 and 34 may be obtained in any suitable manner. 

[0014] The location of any one of coils 32, 33 and 34 can be 
determined in standard ways if the amplitudes of the magnetic signals 
produced by the coil at a sufficient number of detectors 24 is known. 
10 This invention provides a novel noise-insensitive mechanism for finding 
these amplitudes. The detectors 24 used for locating the coils may be the 
same detectors used to obtain MEG data. The detectors may be used 
simultaneously to obtain MEG data and to obtain information used to 
locate the coils. 

15 

[0015] Detectors 24 are sampled periodically to provide a sequence 
of samples. The sequence of samples is acquired at a sampling rate 
sufficient to detect the signals from the coils. The sampling frequency is 
typically at least twice, and preferably at least four times the highest coil 

20 frequency. Sampling is performed during an integration time, T. The 
integration time Tis selected so the signals from the coils 32, 33, 34 can 
be separated from each other and from background noise. Ideally, when 
the output of detectors 24 is integrated over the integration time, signals 
at the noise frequency will integrate to zero over the integration time. 

25 The integration time is selected based upon the expected frequency of 
noise. 
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[0016] Since the position of a subject's head must be updated 
frequently, the integration time must be short, and this results in "spectral 
leakage" where the measurement of a head-coil signal at one frequency is 
affected by signals from the other head coils at other frequencies and/or 
5 by background noise. Spectral leakage can be reduced by multiplying 
the outputs of detectors 24 by a taper function in the integrations, and by 
selecting coil frequencies to be frequencies that do not interfere with 
each other. 

10 [0017] Since the detectors are sampled over finite integration times, 
the locations of the zeros in the noise power spectrum depends upon the 
shape of the taper function, if any, applied to the samples. The examples 
discussed herein relate to the special case where the taper function is a 
"boxcar" (i.e. all of the samples in the integration time Tare weighted 

15 equally) and the noise environment is dominated by signals at a single 
frequency / NO ise an ^ * ts harmonics. The locations of the zeros will differ 
if other taper functions are used. If a taper function other than a boxcar is 
used then the coil frequencies may be selected to be at zeros of the power 
spectrum of the expected noise signal convolved with the taper function. 

20 The issue of separating signals and spectral leakage with finite 

integration times T is discussed in many signal-analysis references, for 
example "Random data: Analysis and Measurement Procedures" by J.S. 
Bendat and A.G.Piersol. 

25 [0018] Figure 2 is a flowchart illustrating a method 100 according 
to one embodiment of the invention. After method 100 begins at block 
102, magnetic signals are generated at block 104. The magnetic signals 
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are generated by means of magnetic emitters, which preferably comprise 
coils. At block 106, the magnetic signals and the noise are measured by 
detectors 24. At block 108, the relative amplitudes of the magnetic 
signals are determined, as described further below. 

5 

[001 9] Where the frequency of the noise signal is f NO isE> an ^ a 
boxcar taper function is used, then suitable integration times T satisfy the 
formula: 

N\ 

T=~ (1) 

/ NOISE 

10 where Nj is an integer and N 2 >2. For example, if the noise fundamental 
frequency / NO ise is 60 Hz then the integration time may be selected to be 
TV/60 seconds. If N 2 =4 then an integration time of 1/15 second is used. If 
the noise fundamental frequency is 50 Hz then the integration time may 
be selected to be TV/50 seconds. If N;=5 an integration time of 1/10 

15 second is used. The integration time is preferably chosen to be short 
relative to any likely motions of the subject's head. 



[0020] The coil frequencies f h f 2 and f 3 are selected to be different 
from one another and at or near frequencies to which the noise 
20 fundamental frequency and its harmonics cannot spread, which are 

referred to herein as ideal coil frequencies f ql . Ideally, each coil frequency 
f qi satisfies the formula: 

f = Ni = Ni x f NOISE 
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Where N 2 is an integer that is not a multiple of 7V 7 . For example, where 
the noise fundamental frequency is 60 Hz and Nj=2 then Equation (2) 
yields frequencies 30Hz, 90Hz, 150Hz, 210Hz ... etc. In preferred 
embodiments, the coil frequencies are chosen so that they will not 
5 interfere with MEG measurements being made. It is also desirable to 
maintain the coil frequencies low enough that currents induced in the 
subject by the fluctuating magnetic fields of the coils are too small to 
interfere significantly with the coil signals. For example, in some MEG 
studies, the MEG signals of interest have frequencies of about 10Hz. In 

10 such cases, the coil frequencies may be chosen to be sufficiently large 
that the lowest coil frequency is significantly higher than the frequencies 
of MEG signals to be studied. Typically the coil frequencies are in the 
range of 25 Hz to 275 Hz. Higher coil frequencies could also be used. 
As another example, if a user is interested in MEG signals having 

15 frequencies in the range of 70 -330Hz, the user could choose Nj =8 and 
N 2 = 3, 5, 7. This choice of N 2 yields coil frequencies of 22.5, 37.5 and 
52.5Hz. In this example, the user would need to filter power line noise 
signals (including harmonics) which are within the range of frequencies 
of the MEG signals being studied. 



20 



[0021] The complex amplitude, A q m , of the signal from the q th coil 
at the m th detector 24 can be determined from the set of samples S k m 
acquired over an integration period at the m th detector by computing the 
sum: 



25 A™ = — 
q N 



cos(2^ / /,)- sin(2^ / f s ) 

L*=o *=o 



= R nm + il m (3) 
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where/^ is the frequency of the 0 th coil, f s is the detector sample rate, 
N=Tf s and k is an index indicating the order of the samples in the 
sequence of samples from each of the magnetic detectors. For locating 
the coils it is only necessary to obtain the relative amplitudes (including 
5 sign) of the signals from each coil at the different detectors (the complex 
phase of the signal from each coil is not important). A simple way to 
obtain the relative amplitudes of the signals is to use a detector which 
detects a strong signal (for example, the detector which detects the 
strongest signal) to establish the phase. If detector m 0 has the strongest 
10 signal at frequency f q then the real relative amplitudes of the signals can 
be given by: 



R; = Re 



AT x 



4 m0 



q " A m0 

A * ) 



(4) 



[0022] Equation (4) is sufficiently accurate as long as the coil 
signals detected by at least a few of detectors 24 are strong enough that 
15 noise does not cause significant errors in the phases of the detected 
signals. More sophisticated techniques for obtaining the relative real 
amplitudes of the coil signals could be used but are usually unnecessary. 
For example, a least squares best estimate of the phase, <j> , and relative 

real amplitude, R q m , of the coil signal of coil q could be obtained as 
20 follows: 

lY IR 

m 



_ qm 

tan K) = W^ 7T\ (S) 
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and, 

R; = Re((^ m + il qm y+) (6) 

[0023] In some cases one or more of the frequencies f h f 2 and f 3 of 
the coil signals may not be exactly equal to the ideal frequencies f qi 
5 identified by Equation (2). This could occur, for example, if the circuits 
used to drive the coils provide limited choice of coil frequencies. In such 
cases one should still use the ideal frequencies f qi in Equation (3) and not 
the actual frequencies of the coil signals. To obtain the best accuracy in 
cases where the frequencies f h f 2 and f 3 of the coil signals are not exactly 
10 equal to the ideal frequencies identified by Equation (2) one should 
compensate for spectral leakage from the signals of the other coils. 

[0024] In cases where one or more of the coil frequencies is not 
ideal, it is desirable that the actual frequencies of the coil signals be 
1 5 related to ideal frequencies of the coil signals by the following 
inequality: 

- A) 2 + U* ~ A) 2 + A ~ A) 2 * if ( ? > 

20 Where equation (7) is satisfied, the actual coil frequencies are 

"substantially equal" to the ideal coil frequencies. It can be seen that 
equation (7) will be satisfied if none of the coil frequencies departs from 

03 

a corresponding ideal frequency by more than approximately — . 
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[0025] To compensate for spectral leakage that occurs when non- 
ideal coil frequencies are used, one can define the functions G c and G s as 
follows: 

T 

G c {f q ,f qi ) = | i\cos{2nf q t)e n ^ dt 

1 2 

sin«/, - f„)T) | sin(4/, + f ¥ )T) (8) 

*(/,-/„)*■ + 4f< + f«) T 



and, 



= f J\sin(2^r)e-' 2 ^Vr 

2 

(9) 



10 



where f qi is the ideal frequency closest to f r The real and imaginary parts 
C q and S q of the amplitudes measured at ideal frequencies f qi can be 
related to the real and imaginary parts R q and I q of the true head coil 
amplitudes by way of the matrix equations: 



G c {A,fu) G c (f 2 ,f u ) G c (f 3 ,f u ) 
G c {A,f 2i ) G c (f 2 ,f 2i ) G c (f 3 ,f 2i ) 
Gclfi.f*) G c (f 2 ,f 3i ) G c (f 3 ,f 3 ) 







X 









(10) 



and, 



A. 



G s {A,/ u ) G s {f 2 ,f u ) G s (f 3 ,f u ) 
G s (A,f 2i ) G s (f 2 ,f 2 ) G s (f 3 ,f 2i ) 
G s {A,f v ) G s (f 2 ,f 3i ) G s (f 3 ,f 3i ) 



(11) 
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The simple form of equations (10) and (1 1) as a pair of 3x3 matrix 
equations can be achieved by expressing the integration time as the 
interval -772 < t < T/2. Equations (10) and (11) can be replaced with an 
equivalent single 6><6 matrix equation or an equivalent single 3x3 matrix 
5 equation in which the variables and coefficients have complex values. 
Equations (10) and (1 1), or their mathematical equivalents, may be 
applied to obtain corrected amplitudes to compensate for the coil 
frequencies not being exactly at the ideal frequencies described above. 

[0026] Motion of the head coils will result in a time-varying signal 
amplitude, which in turn will cause spectral leakage such that the signal 
from one head coil affects the measurement of the amplitude of the 
signals from the other head coils, even if the coils are operated at the 
ideal frequencies specified in equation (2). Where the integration time is 
short enough that the signal amplitude can be modelled as a quadratic 
function, the spectral leakage can be considered to be a function of 
frequency separation only. Spectral leakage due to motion of the coils 
decreases with increased separation between the frequencies f h f 2 and f 3 . 
It is desirable to maintain a spacing of, for example, 30 Hz or more, 
between all pairs of f h f 2 and f 3 to minimize spectral leakage. 

[0027] It can be appreciated that selecting coil frequencies as 
described above is a special case of a more general method. In general, 
the coil frequencies are selected to be equal to or substantially equal to 
25 those frequencies at which the power spectrum of the noise signal is zero. 
In this context, "substantially equal" is defined in equation (7) above. 



10 



15 



20 
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0.3 

Frequencies which are all within about ~jr , where Tis the integration 

time, of the corresponding ideal frequencies are substantially equal to the 
ideal frequencies. This is typically close enough to compensate for 
errors which arise from the coil signals not being at the ideal frequencies 
5 by using, for example, equations (10) and (11). 

[0028] Standard techniques may be used to determine the locations 
of each of the coils based upon the amplitudes determined above. Some 
suitable head localization algorithms are described in de Munck et al., 

10 Phys. Med. Biol. 46 (2001) pp. 2041-2052. For example, one way to 
determine the positions of the coils can be used in the case in which the 
coils are each small enough to be considered a point source and the 
detectors are magnetometers or first order gradiometers. It can be shown 
that the position y can be determined by finding a value for y which 

1 5 minimizes the function H where: 

H(y)= Z (^)'QW" s k ) 2 = mAfh- 2b -m+ £ s 2 k (12 ) 

k k 

where 

^00 is the best-fit estimate of the source magnetic moment if the 
20 source is at point y . It is related to matrix A and vector b by 

fh[y) = A~ x b 

Matrix A and vector b are functions of y and are defined below; 
s k is the amplitude of a the signal from a coil in the k th detector; 
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Ck M (y ) is th e response of sensor A: to a magnetic source located at 
point y with moment of unit strength and oriented in spatial 
direction ja. C k (y) may be represented as a power series in y 

and the parameters that describe sensor k, or it may be expressed 
5 as an analytic function of the same variables; 

Matrix A is a function of source point y . The components of A 

k 

Vector b is a function of source point y . The components of b are 
K = Z s k c k»{y)- 

k 

10 

[0029] It is not necessary to use the signals detected at all detectors 
to determine the locations of each of the coils. It is generally acceptable 
to use signals detected at a subset of the detectors. Acceptable results can 
generally be obtained using coil signals detected at 30 or fewer properly 

1 5 chosen detectors. In some embodiments of the invention coil signals 

from 12 to 15 detectors are used to determine the position of each coil. It 
is not necessary that the same subset of coils be used to determine the 
positions of all of the coils. Indeed, it is desirable to use a different 
subset of the detectors for determining the position of each of the coils. 

20 Since the coils are expected to be in positions which are roughly known 
one could establish a predetermined subset of detectors to use for each 
coil. 
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[0030] A subset of the magnetic detectors to use for determining 
the position of a coil may be identified in various ways. In general it is 
desirable to include in the subset detectors having outputs which vary 
strongly with the position of the coil relative to the detector. In some 
5 embodiments of the invention multiple criteria are used to select 
detectors in the subset. For example, in some embodiments, some 
detectors are included in the subset because they have relatively large 
signals from the coil and other detectors are included in the subset 
because the coil signal that they detect has a large variation with coil 
1 0 position. 



[0031] One algorithm for selecting detectors to include in the 
subset for a coil involves the following steps: 

Determine a number M f of detectors to use for the fit. M f may be 
1 5 predetermined. M f in the range of 12 to 15 is generally sufficient. 

must be at least 6. More detectors may be required for coils 
which are relatively far from the closest detectors. 
Sort the detectors in order of the amplitude of the signal detected 
from the coil. Select the detectors having the largest signals for 
20 inclusion in the subset. may be, for example, 27V/3. In general 

may be any significant fraction of M f (i.e. a fraction greater 
than about 1/37V). 

Calculate an indicator, V, of the variation in the signal with coil 
position for each detector and include in the subset the MfM^ 
25 detectors having the largest values of V. Various functions can be 

used as the indicator. For example, in some embodiments, Kis 
given by: 
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V = 



\ S m S max\ 



r - r 

m max 



(13) 



where s max is the amplitude of the signal at the detector having the 
largest signal from the coil and f max is the location of the detector 
having the largest signal from the coil. 

5 

[0032] Information about the coil positions may be used in various 
ways. For example, the coil positions may be monitored in real time and 
compared to reference positions. Suitable action may be taken based 
upon the coil positions. For example: 
10 an alarm may be triggered if the coil positions are determined to be 

too far away from the reference positions; 

data may be discarded for those times during which the coil 

positions are determined to be too far away from the reference 

positions; 

1 5 • a sensory (e.g. visual, audible or tactile) feedback indicator may be 
provided to help subjects keep their heads in a desired position; 
• the positions of the coils may be recorded along with MEG data 
for later correction of the MEG data to account for motion of the 
subject's head; 

20 • the MEG data may be corrected in real time to correct for motions 
of the subject's head. One way to correct the MEG data is 
described in the co-pending commonly owned patent application 
entitled MOTION COMPENSATION IN BIOMAGNETIC 
MEASUREMENT being filed together herewith, which is hereby 

25 incorporated by reference herein. 
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[0033] Certain implementations of the invention comprise computer 
processors which execute software instructions which cause the 
processors to perform a method of the invention. For example, one or 
5 more processors in a controller for a MEG system may implement a 
method of Figure 2 by executing software instructions from a program 
memory accessible to the processor(s). The invention may also be 
provided in the form of a program product. The program product may 
comprise any medium which carries a set of computer-readable signals 

10 comprising instructions which, when executed by a computer processor, 
cause the data processor to execute a method of the invention. Program 
products according to the invention may be in any of a wide variety of 
forms. The program product may comprise, for example, physical media 
such as magnetic data storage media including floppy diskettes, hard disk 

15 drives, optical data storage media including CD ROMs, DVDs, electronic 
data storage media including ROMs, flash RAM, or the like or 
transmission-type media such as digital or analog communication links. 

[0034] Where a component (e.g. a software module, processor, 
20 detector, assembly, device, circuit, etc.) is referred to above, unless 

otherwise indicated, reference to that component (including a reference 
to a "means 1 ') should be interpreted as including as equivalents of that 
component any component which performs the function of the described 
component (i.e., that is functionally equivalent), including components 
25 which are not structurally equivalent to the disclosed structure which 
performs the function in the illustrated exemplary embodiments of the 
invention. 
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[0035] As will be apparent to those skilled in the art in the light of 
the foregoing disclosure, many alterations and modifications are possible 
in the practice of this invention without departing from the spirit or scope 
thereof For example: 

• In the above description, a number of coil signals are monitored 
simultaneously. In alternative embodiments of the invention the 
coil signals are time multiplexed. In this case, the coil signals may 
be at the same frequency or at different frequencies. The coil 
frequency or frequencies are chosen to be frequencies at which the 
power spectrum of the expected noise signal convolved with any 
taper function being applied has zeros. 

• The foregoing discussion assumes that the magnetic signals 
generated by the coils are each at a single fixed frequency and 
remain uniform in frequency and amplitude over the integration 
period. In general, this is not necessary. 

Magnetic detectors 24 may be of any suitable types. For example, 
detectors 24 could be magnetometers or first or second order 
magnetic gradiometers. Detectors 24 could be SQUID 
magnetometers or gradiometers. 

• While the invention is discussed above in the context of a MEG 
system, the invention may be applied to locating other objects 
oscillating magnetic fields relative to an array of magnetic sensors. 

• While the mathematical formulae presented herein have been 
structured to avoid the need for computations using complex 
numbers, one could perform equivalent calculations using complex 
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numbers at the expense of somewhat greater memory usage and 

somewhat slower calculation. 
Accordingly, the scope of the invention is to be construed in accordance 
with the substance defined by the following claims. 



